In [1]:
#import functions
%pylab inline
# from MyUnits import *
from MyFunctions import *
from qutip import *
# from MyQubit import *
# import mpld3
import multiprocessing as mp
import itertools
In [2]:
import scipy.constants as sc
In [3]:
import time
In [4]:
def calc_spectrum_3(N, M, Delta_C,Delta_NR,
Delta_Q, g, lamb,
A_C, A_Q, A_NR,
kappa_n, gamma_rel,
gamma_dep,n_th_a):
"""
Calculate Qubit density Matrix rho and trace of: rho_sz, rho_n, rho_x, rho_a : Based on work of : Suri, B. et al. Observation of Autler–Townes effect in a dispersively dressed Jaynes–Cummings system. New J. Phys. 15, 125007 (2013).
look also calc_spectrum function
N: number states cavity
t_list: time - not used
Delta_d_tilde: cavity tone frequency difference
Delta_s_tilde: qubit tone frequency difference
chi: dispersive coupling
Gamma_d: cavity tone power (GHz)
Gamma_s: qubit tone power (GHz)
"""
# cavity operators
a = tensor(destroy(N), qeye(M), qeye(2))
n_C = a.dag() * a
x_C = a + a.dag()
# NR operators
b = tensor(qeye(N), destroy(M), qeye(2))
n_NR = b.dag() * b
x_NR = b + b.dag()
# qubit operators
sm = tensor(qeye(N), qeye(M), destroy(2))
sz = tensor(qeye(N), qeye(M), sigmaz())
sx = tensor(qeye(N), qeye(M), sigmax())
nq = sm.dag() * sm
xq = sm + sm.dag()
# Unitary
I = tensor(qeye(N), qeye(M), qeye(2))
# Hamiltonian
H_C = Delta_C * (a.dag() * a + I / 2.0) # Cavity
H_NR = Delta_NR * (b.dag() *b + I / 2.0 ) # Nanoresonator
H_Q = (Delta_Q / 2.0) * sz # qubit
H_QtoC = g * (sm * a.dag() + sm.dag()*a ) # interaction Cavity Qubit
H_QtoNR = lamb * (b.dag() + b)*sx # interaction NR Qubit
H_A = A_C * (a.dag() + a) # tone cavity
H_B = A_Q * (sm.dag() + sm) # tone qubit
H_C = A_NR * (b.dag() + b) # tone NR
# # extra term Unitary tranformation and RWA
# He = nc * chi ** 2 / (4 * Delta) * (w_q / Delta - I)
# Hd = (Gamma_d / 2.0) * (a + a.dag()) # Drive cavity
# Hs = (Gamma_s / 2.0) * (sm + sm.dag()) # Drive qubit
# HKerr = zeta_1 * ((a.dag() * a) ** 2) * sz + zeta_2 * (a.dag() * a) ** 2 # non linear Kerr terms
H = H_C + H_NR + H_Q + H_QtoC + H_QtoNR + H_A + H_B + H_C # total Hamiltonian
# collapse operators
c_op_list = []
rate = kappa_n * (1 + n_th_a)
if rate > 0.0:
c_op_list.append(sqrt(rate) * a)
rate = kappa_n * n_th_a
if rate > 0.0:
c_op_list.append(sqrt(rate) * a.dag())
rate = gamma_rel * (1 + n_th_a)
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm)
rate = gamma_rel * (n_th_a)
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm.dag())
rate = (gamma_dep / 2) * (1 + n_th_a)
if rate > 0.0:
c_op_list.append(sqrt(rate) * sz)
# Calculation stedy state Master equation
rho = steadystate(H, c_op_list)
# variables to be returned
rho_sz = rho * sz
# rho_n = rho * (a.dag() * a + I / 2)
# rho_x = rho * (a.dag() + a)
rho_a = rho * a
rho_b = rho * b
return rho_sz.tr(), rho_a.tr(), rho_b.tr()
In [5]:
N = 5
M = 5
Delta_C = 1
Delta_NR = 1
Delta_Q = 1
g = 0.119
lamb = 0.1
A_C = 0.000308
A_Q = 0.00026
A_NR = 0.000001
kappa_n = 0.00025
gamma_rel = 0.00023
gamma_dep = 0.00114
n_th_a = 0.3
In [6]:
def Delta(f_r,f_d):
"""
f_r resonance freq
f_d drive freq
"""
return f_r-f_d
def Gamma(P,Q_c,Q_l,kappa,w):
return sqrt(Q_c*kappa*1e9*dBmtoW(P)/(2*Q_l*sc.hbar*w*1e9))/1e9
In [ ]:
In [7]:
# a1, b1 = zip(*itertools.product(Delta_C, Delta_Q))
# pool = mp.Pool(processes=12)
# results2 = [pool.apply_async(calc_spectrum_3,
# (N,
# M,
# C,
# Delta_NR,
# Q,
# g,
# lamb,
# A_C,
# A_Q,
# A_NR,
# kappa_n,
# gamma_rel,
# gamma_dep,
# n_th_a)) for C,Q in zip(a1,b1)]
# results = [p.get() for p in results2]
# results.sort()
In [8]:
Ej = 11.55 # Maximum Josephson Energy
Ec = 0.22 # Capacitive Energy
w_ge_max = sqrt(8*Ec*Ej)-Ec
print('Qubit max frequency: ',w_ge_max, 'GHz')
f = w_ge_max *1e9 # GHz
T = 50e-3 # K
n_th_a = 1/(exp(sc.h*f/(sc.k*T)-1))
print('<n> thermal: ',n_th_a)
In [10]:
Q_c = 20000
Q_l = Q_c
N = 5
M = 5
g = 0.1
lamb = 0.00011
kappa_n = 0.0002
gamma_rel = 0.0002
gamma_dep = 0.001
# Cavity
w_c = 5 * (2 * pi)
w_d = w_c + 0.006283185307175643 #(2 * pi) * 4.995
P_C = -120# linspace(-60,-150,20)
# Nanoresonator
w_NR = 3 * (2 * pi)
w_p = (2 * pi) * 2.9996
P_NR = -1000
# Qubit
w_q = 4.5 * (2 * pi)
w_s = (2 * pi) * linspace(4.2,4.8,300)#4.501
P_Q = linspace(-50,-150,50)
In [11]:
Delta_C = Delta(w_c, w_d)
Delta_NR = Delta(w_NR, w_p)
Delta_Q = Delta(w_q, w_s)
# print(Delta_C,Delta_NR,Delta_Q)
In [12]:
A_C = Gamma(P_C,Q_c,Q_l,kappa_n,w_c)
A_Q = Gamma(P_Q,Q_c,Q_l,kappa_n,w_c)
A_NR = Gamma(P_NR,Q_c,Q_l,kappa_n,w_c)
# print(r'A_C = ', A_C, '\n\nA_Q = ', A_Q, '\n\nA_NR = ' , A_NR)
In [13]:
a1, b1 = zip(*itertools.product(A_Q, Delta_Q))
In [ ]:
pool = mp.Pool(processes=12)
t0 = time.time()
results1 = [pool.apply_async(calc_spectrum_3,
(N,
M,
Delta_C,
Delta_NR,
y,
g,
lamb,
A_C,
x,
A_NR,
kappa_n,
gamma_rel,
gamma_dep,
n_th_a)) for x , y in zip(a1,b1)]
t1 = time.time()
print("elapsed =", (t1-t0))
# results2 = [pool.apply_async(calc_spectrum_3,
# (N,
# M,
# Delta_C,
# Delta_NR,
# Delta_Q,
# g,
# lamb,
# A_C,
# A_Q
# A_NR,
# kappa_n,
# gamma_rel,
# gamma_dep,
# n_th_a)) for x , y in zip(a1,b1)]
results = [p.get() for p in results1]
In [ ]:
results = [p.get() for p in results1]
In [ ]:
shape(results)
In [ ]:
Tr_sz = reshape(array(results)[:,0],(-1,len(Delta_Q+1)))
Tr_a = reshape(array(results)[:,1],(-1,len(Delta_Q+1)))
Tr_b = reshape(array(results)[:,2],(-1,len(Delta_Q+1)))
In [ ]:
fig, ax = subplots(2,2, figsize=(16,10))
im = ax[0,0].pcolor( w_s/2/pi,P_Q,abs(Tr_sz))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[0,0])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[0,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[0,0].set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)
im2 = ax[0,1].pcolor(w_s/2/pi,P_Q,abs(Tr_a))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[0,1])
# ax[0,1].set_ylim(P_i,P_f)
# ax[0,1].set_xlim(4.27,4.39)
ax[0,1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[0,1].set_title(r'$Tr[\rho {a}]$',fontsize=20)
im3 = ax[1,0].pcolor(w_s/2/pi,P_Q,abs(Tr_b))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1,0])
# ax[1,0].set_ylim(P_i,P_f)
# ax[1,0].set_xlim(4.27,4.39)
ax[1,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[1,0].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[1,0].set_title(r'$Tr[\rho (b)]$',fontsize=20)
In [ ]: